helper function to process time series

And some package and labeling info

library(ggplot2)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)
library(factoextra) # easier visualizing outputs of prcomp
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
timeseries_dir <- 'extracted_timeseries/'

metrics_write_dir <-  'extracted_timeseries/extracted_metrics/'


# get ecs ordering/labels
esm_labels <- read.csv(paste0(timeseries_dir,'global_tas_allesms.csv'), stringsAsFactors = FALSE) %>%
  select(esm) %>% distinct %>% 
  mutate(plotesm = paste0(letters[as.integer(row.names(.))], '.', esm),
         ECS_order = as.integer(row.names(.)))
print(esm_labels)
##              esm         plotesm ECS_order
## 1    CAMS-CSM1-0   a.CAMS-CSM1-0         1
## 2         MIROC6        b.MIROC6         2
## 3      GFDL-ESM4     c.GFDL-ESM4         3
## 4      FGOALS-g3     d.FGOALS-g3         4
## 5  MPI-ESM1-2-HR e.MPI-ESM1-2-HR         5
## 6  MPI-ESM1-2-LR f.MPI-ESM1-2-LR         6
## 7     MRI-ESM2-0    g.MRI-ESM2-0         7
## 8  ACCESS-ESM1-5 h.ACCESS-ESM1-5         8
## 9   IPSL-CM6A-LR  i.IPSL-CM6A-LR         9
## 10   CESM2-WACCM   j.CESM2-WACCM        10
## 11   UKESM1-0-LL   k.UKESM1-0-LL        11
## 12       CanESM5       l.CanESM5        12

Spatial info

we want the shapes for plotting anyway, so prep them

## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
## Reading layer `IPCC-WGI-reference-regions-v4' from data source 
##   `/Users/snyd535/Documents/GitHub/stitches_in_r/R/inst/shinyApp/python_curation/IPCC-WGI-reference-regions-v4_shapefile/IPCC-WGI-reference-regions-v4.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 58 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: -90 xmax: 180 ymax: 90
## Geodetic CRS:  WGS 84

Load ESM data that’s been processed

region_summary_main <- read.csv(paste0(metrics_write_dir, 'IPCC_land_regions_metrics.csv'), 
                                stringsAsFactors = FALSE)
region_summary_main %>%
  filter(experiment != 'ssp119',
         experiment != 'ssp434',
         experiment != 'ssp460',
         experiment != 'ssp534-over')  %>%
    rename(region = acronym) ->
  region_summary 

print(head(region_summary))
##          esm experiment ensemble variable       type              name region
## 1 ACCESS-CM2     ssp126 r1i1p1f1       pr       Land Arabian-Peninsula    ARP
## 2 ACCESS-CM2     ssp126 r1i1p1f1       pr       Land    Central-Africa    CAF
## 3 ACCESS-CM2     ssp126 r1i1p1f1       pr Land-Ocean         Caribbean    CAR
## 4 ACCESS-CM2     ssp126 r1i1p1f1       pr       Land       C.Australia    CAU
## 5 ACCESS-CM2     ssp126 r1i1p1f1       pr       Land   C.North-America    CNA
## 6 ACCESS-CM2     ssp126 r1i1p1f1       pr       Land            E.Asia    EAS
##           iasd      end_val   end_anomaly end_anomaly_pct mid_century_val
## 1 7.227917e-07 1.700625e-06  4.306991e-07     0.339153029    2.051770e-06
## 2 3.009324e-06 4.828026e-05 -4.965491e-08    -0.001027416    4.941042e-05
## 3 6.106140e-06 3.052769e-05  2.089845e-07     0.006892922    3.429079e-05
## 4 2.830779e-06 1.130480e-05 -1.269868e-06    -0.100986213    1.126320e-05
## 5 2.861294e-06 3.232150e-05  1.953880e-06     0.064340892    3.222543e-05
## 6 2.687448e-06 5.449186e-05  5.826705e-06     0.119730531    5.260652e-05
##     mid_anomaly mid_anomaly_pct
## 1  7.818446e-07      0.61566180
## 2  1.080502e-06      0.02235680
## 3  3.972086e-06      0.13101110
## 4 -1.311461e-06     -0.10429388
## 5  1.857811e-06      0.06117738
## 6  3.941364e-06      0.08098944

Prep data for PCA

  • idea is that there’s something really similar across ESMs, especially for SSP126/245 -> their RGB maps have so much in common.

  • is there a sort of core map that gets deviations added to as you move through scenarios and ESMs? Or how many distinct maps do we actually have

  • And the aim of the archive is to establish that our selection of ESMs forms a spanning set of vectors for all/most of CMIP6

  • observation = scenarioXesm

  • individual variable is a regionXmetric

  • so each observation vector is 3*Nregions long: <Tr1...Trn, Pr1...Prn, IASDr1...IASDrn > for a scenarioXesm <-> there are 48 observations for each of 129 variables

  • And then becaues we’re interested in a pseudo-basis for the CMIP6 archive, we want the training data to be all scenarios and ESMs

  • Out of sample tests -> how does basis do on overshoots that it doesn’t train on; withhold one of the 12 ESMs and learn from other 11 to see how it goes; then see how a completely external ESM does on training from all 12

  • nothing a priori to preclude us from doing EOFs on the gridded time series data, or the gridded metrics, but having on IPCC regions avoids the issue of all ESMs being on their own grids

  • and overall doing region metrics helps keep data size operating on manageable

# reshape:
grouped_data <- split(region_summary, list(region_summary$esm, region_summary$experiment))
grouped_data <- grouped_data[sapply(grouped_data, function(x) dim(x)[1]) > 0]

shaped_data <- lapply(grouped_data, function(group){
    if (nrow(group) >0) {
      group %>%
        filter(variable == 'tas') %>% 
        select(esm, experiment, ensemble,type, region,
               iasd, end_anomaly, mid_anomaly) %>%
        rename(tas_iasd = iasd,
               tas_end_anomaly = end_anomaly,     
               tas_mid_anomaly = mid_anomaly) %>% 
        left_join(group %>%
                      filter(variable == 'pr') %>% 
                      select(esm, experiment, ensemble,type, region,
                             iasd, end_anomaly_pct, mid_anomaly_pct) %>%
                      rename(pr_iasd = iasd,
                             pr_end_anomaly_pct = end_anomaly_pct,    
                             pr_mid_anomaly_pct = mid_anomaly_pct),
                  by = c('esm','experiment', 'ensemble', 'type', 'region')
               ) ->
    reshaped
    
    reshaped %>%
        group_by(esm, experiment, type, region) %>%
        summarize(tas_iasd=mean(tas_iasd),
                  tas_end_anomaly=mean(tas_end_anomaly, na.rm = T),
                  tas_mid_anomaly = mean(tas_mid_anomaly, na.rm = T),
                  pr_iasd = mean(pr_iasd, na.rm = T),
                  pr_end_anomaly_pct = mean(pr_end_anomaly_pct, na.rm = T),
                  pr_mid_anomaly_pct = mean(pr_mid_anomaly_pct, na.rm = T) ) %>%
        ungroup() %>%
        mutate(ensemble = 'ensemble_avg') -> #%>%
        #bind_rows(reshaped) ->
        reshaped2
    
    #TODO need to change shaping if including ensemble members
    reshaped2 %>%
        select(-type) %>%
        gather(metric, value, -esm, -experiment, -ensemble, -region) %>%
        mutate(row_id = paste0(region, '~', metric),
           col_id = paste0(esm, '~', experiment, '~', ensemble)) %>%
        select(-esm, -experiment,-ensemble, -region, -metric) %>%
        as.data.frame() ->
    reshaped3 
  
  colnames(reshaped3) <- c(paste0(unique(reshaped3$col_id)), 'row_id', 'col_id')
  rownames(reshaped3) <- paste0(reshaped3$row_id)
  
  reshaped3 %>% 
    select(-col_id, -row_id) ->
    out
  
  return(out)  
    }
    
    }
)
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
# combine columns but then transpose because prcomp wants rows to be observations 
# and columns to be variables (but doing it the other way first is easier to code)
data <- t(do.call(cbind, shaped_data) )

print(str(data))
##  num [1:111, 1:258] 0.331 0.291 0.278 0.323 0.311 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : chr [1:111] "ACCESS-CM2~ssp126~ensemble_avg" "ACCESS-ESM1-5~ssp126~ensemble_avg" "BCC-CSM2-MR~ssp126~ensemble_avg" "CAMS-CSM1-0~ssp126~ensemble_avg" ...
##   ..$ : chr [1:258] "ARP~tas_iasd" "CAF~tas_iasd" "CAU~tas_iasd" "CNA~tas_iasd" ...
## NULL

do PCA

have to remove some NA’s:

  • manually checked all 4 - no PR ens average values
data1 <- na.omit(data)
print('removed rows due to missing data')
## [1] "removed rows due to missing data"
print(row.names(data)[c(which(!(row.names(data) %in% row.names(data1))))])
## [1] "CIESM~ssp126~ensemble_avg"      "CIESM~ssp245~ensemble_avg"     
## [3] "CIESM~ssp585~ensemble_avg"      "NorESM2-LM~ssp585~ensemble_avg"
data_pca <- prcomp(data1, center=TRUE, scale = TRUE)

# quick visualize
factoextra::fviz_eig(data_pca)

further visualization

# coordinates of each ESM-SSP combo in the PCs
data_pca$x %>%
  as.data.frame() %>%
  mutate(id = row.names(.)) %>%
  separate(id, into=c('model', 'scenario', 'ensemble'), sep = '~') %>%
  select(model, scenario, PC1, PC2, PC3, PC4, PC5, PC6, PC7, PC8, PC9)->
  coordinates

PC1 vs PC2

ggplot(coordinates) +
    geom_point(mapping = aes(x = PC1, y = PC2, color = model, shape = scenario)) 

ggplot(coordinates) +
    geom_text(mapping = aes(x = PC1, y = PC2, color = model, label = model), size=0.95) 

ggplot(coordinates) +
  geom_point(mapping = aes(x = PC1, y = PC2, color = scenario))

PC1 vs PC3

ggplot(coordinates) +
  geom_point(mapping = aes(x = PC1, y = PC3, color = model, shape = scenario))

ggplot(coordinates) +
  geom_point(mapping = aes(x = PC1, y = PC3, color = scenario))

PC1 vs PC4

ggplot(coordinates) +
  geom_point(mapping = aes(x = PC1, y = PC4, color = model, shape = scenario))

ggplot(coordinates) +
  geom_point(mapping = aes(x = PC1, y = PC4, color = scenario))

PC2 vs PC3

ggplot(coordinates) +
  geom_point(mapping = aes(x = PC2, y = PC3, color = model, shape = scenario))

ggplot(coordinates) +
  geom_point(mapping = aes(x = PC2, y = PC3, color = scenario))

PC2 vs PC4

ggplot(coordinates) +
  geom_point(mapping = aes(x = PC2, y = PC4, color = model, shape = scenario))

ggplot(coordinates) +
  geom_point(mapping = aes(x = PC2, y = PC4, color = scenario))

PC3 vs PC4

ggplot(coordinates) +
  geom_point(mapping = aes(x = PC3, y = PC4, color = model, shape = scenario))

ggplot(coordinates) +
  geom_point(mapping = aes(x = PC3, y = PC4, color = scenario))

plot PCs as maps

pcs <- data_pca$rotation

str(pcs) 
##  num [1:258, 1:107] 0.0179 0.0129 0.0288 0.0218 0.0151 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : chr [1:258] "ARP~tas_iasd" "CAF~tas_iasd" "CAU~tas_iasd" "CNA~tas_iasd" ...
##   ..$ : chr [1:107] "PC1" "PC2" "PC3" "PC4" ...
# pull the region and metric info back in
as.data.frame(pcs) %>%
  mutate(row_id = row.names(.)) %>%
  separate(row_id, into = c('region', 'metric'), sep ='~') ->
  tmp


# We want to convert each set of 3 metrics in a region to a single
# rgb coded value.
# and we want all on same color scale. 

# so reshape for each metric (which goes to one of r, g, b)
# with all PCs want to consider stacked.
tmp %>%
  filter(metric == 'tas_end_anomaly') %>%
  select(region, metric, PC1) %>%
  rename(tas_end_anomaly = PC1) %>%
  mutate(pc = 'PC1') %>%
    bind_rows(tmp %>%
                filter(metric == 'tas_end_anomaly') %>%
                select(region, metric, PC2) %>%
                rename(tas_end_anomaly = PC2) %>% 
                mutate(pc = 'PC2'), 
              tmp %>%
                filter(metric == 'tas_end_anomaly') %>%
                select(region, metric, PC3) %>%
                rename(tas_end_anomaly = PC3) %>% 
                mutate(pc = 'PC3') ,
              tmp %>%
                filter(metric == 'tas_end_anomaly') %>%
                select(region, metric, PC4) %>%
                rename(tas_end_anomaly = PC4) %>% 
                mutate(pc = 'PC4'),
              tmp %>%
                filter(metric == 'tas_end_anomaly') %>%
                select(region, metric, PC5) %>%
                rename(tas_end_anomaly = PC5) %>% 
                mutate(pc = 'PC5'),
              tmp %>%
                filter(metric == 'tas_end_anomaly') %>%
                select(region, metric, PC6) %>%
                rename(tas_end_anomaly = PC6) %>% 
                mutate(pc = 'PC6'),
              tmp %>%
                filter(metric == 'tas_end_anomaly') %>%
                select(region, metric, PC7) %>%
                rename(tas_end_anomaly = PC7) %>% 
                mutate(pc = 'PC7'),
              tmp %>%
                filter(metric == 'tas_end_anomaly') %>%
                select(region, metric, PC8) %>%
                rename(tas_end_anomaly = PC8) %>% 
                mutate(pc = 'PC8'),
              tmp %>%
                filter(metric == 'tas_end_anomaly') %>%
                select(region, metric, PC9) %>%
                rename(tas_end_anomaly = PC9) %>% 
                mutate(pc = 'PC9'),
              tmp %>%
                filter(metric == 'tas_end_anomaly') %>%
                select(region, metric, PC10) %>%
                rename(tas_end_anomaly = PC10) %>% 
                mutate(pc = 'PC10')
              ) ->
  tas
print(head(tas))
##                         region          metric tas_end_anomaly  pc
## ARP~tas_end_anomaly...1    ARP tas_end_anomaly      0.09396300 PC1
## CAF~tas_end_anomaly...2    CAF tas_end_anomaly      0.09282450 PC1
## CAU~tas_end_anomaly...3    CAU tas_end_anomaly      0.09292208 PC1
## CNA~tas_end_anomaly...4    CNA tas_end_anomaly      0.09335706 PC1
## EAS~tas_end_anomaly...5    EAS tas_end_anomaly      0.09580033 PC1
## EAU~tas_end_anomaly...6    EAU tas_end_anomaly      0.09210200 PC1
print(tail(tas))
##                            region          metric tas_end_anomaly   pc
## WNA~tas_end_anomaly...425     WNA tas_end_anomaly     0.048192076 PC10
## WSAF~tas_end_anomaly...426   WSAF tas_end_anomaly     0.004906223 PC10
## WSB~tas_end_anomaly...427     WSB tas_end_anomaly     0.047285891 PC10
## CAR~tas_end_anomaly...428     CAR tas_end_anomaly     0.052966268 PC10
## MED~tas_end_anomaly...429     MED tas_end_anomaly     0.054653960 PC10
## SEA~tas_end_anomaly...430     SEA tas_end_anomaly     0.018720644 PC10
tmp %>%
  filter(metric == 'pr_end_anomaly_pct') %>%
  select(region, metric, PC1) %>%
  rename(pr_end_anomaly_pct = PC1) %>%
  mutate(pc = 'PC1') %>%
  bind_rows(tmp %>%
              filter(metric == 'pr_end_anomaly_pct') %>%
              select(region, metric, PC2) %>%
              rename(pr_end_anomaly_pct = PC2) %>%
              mutate(pc = 'PC2'), 
            tmp %>%
              filter(metric == 'pr_end_anomaly_pct') %>%
              select(region, metric, PC3) %>%
              rename(pr_end_anomaly_pct = PC3) %>%
              mutate(pc = 'PC3') ,
            tmp %>%
              filter(metric == 'pr_end_anomaly_pct') %>%
              select(region, metric, PC4) %>%
              rename(pr_end_anomaly_pct = PC4) %>%
              mutate(pc = 'PC4') ,
            tmp %>%
              filter(metric == 'pr_end_anomaly_pct') %>%
              select(region, metric, PC5) %>%
              rename(pr_end_anomaly_pct = PC5) %>%
              mutate(pc = 'PC5') ,
            tmp %>%
              filter(metric == 'pr_end_anomaly_pct') %>%
              select(region, metric, PC6) %>%
              rename(pr_end_anomaly_pct = PC6) %>%
              mutate(pc = 'PC6') ,
            tmp %>%
              filter(metric == 'pr_end_anomaly_pct') %>%
              select(region, metric, PC7) %>%
              rename(pr_end_anomaly_pct = PC7) %>%
              mutate(pc = 'PC7') ,
            tmp %>%
              filter(metric == 'pr_end_anomaly_pct') %>%
              select(region, metric, PC8) %>%
              rename(pr_end_anomaly_pct = PC8) %>%
              mutate(pc = 'PC8') ,
            tmp %>%
              filter(metric == 'pr_end_anomaly_pct') %>%
              select(region, metric, PC9) %>%
              rename(pr_end_anomaly_pct = PC9) %>%
              mutate(pc = 'PC9') ,
            tmp %>%
              filter(metric == 'pr_end_anomaly_pct') %>%
              select(region, metric, PC10) %>%
              rename(pr_end_anomaly_pct = PC10) %>%
              mutate(pc = 'PC10') 
              ) ->
  pr


tmp %>%
  filter(metric == 'tas_iasd') %>%
  select(region, metric, PC1) %>%
  rename(tas_iasd = PC1) %>%
  mutate(pc = 'PC1') %>%
  bind_rows(tmp %>%
              filter(metric == 'tas_iasd') %>%
              select(region, metric, PC2) %>%
              rename(tas_iasd = PC2) %>%
              mutate(pc = 'PC2'), 
            tmp %>%
              filter(metric == 'tas_iasd') %>%
              select(region, metric, PC3) %>%
              rename(tas_iasd = PC3) %>%
              mutate(pc = 'PC3') ,
            tmp %>%
              filter(metric == 'tas_iasd') %>%
              select(region, metric, PC4) %>%
              rename(tas_iasd = PC4) %>%
              mutate(pc = 'PC4') ,
            tmp %>%
              filter(metric == 'tas_iasd') %>%
              select(region, metric, PC5) %>%
              rename(tas_iasd = PC5) %>%
              mutate(pc = 'PC5'), 
            tmp %>%
              filter(metric == 'tas_iasd') %>%
              select(region, metric, PC6) %>%
              rename(tas_iasd = PC6) %>%
              mutate(pc = 'PC6'), 
            tmp %>%
              filter(metric == 'tas_iasd') %>%
              select(region, metric, PC6) %>%
              rename(tas_iasd = PC6) %>%
              mutate(pc = 'PC6'), 
            tmp %>%
              filter(metric == 'tas_iasd') %>%
              select(region, metric, PC7) %>%
              rename(tas_iasd = PC7) %>%
              mutate(pc = 'PC7'), 
            tmp %>%
              filter(metric == 'tas_iasd') %>%
              select(region, metric, PC8) %>%
              rename(tas_iasd = PC8) %>%
              mutate(pc = 'PC8'), 
            tmp %>%
              filter(metric == 'tas_iasd') %>%
              select(region, metric, PC9) %>%
              rename(tas_iasd = PC9) %>%
              mutate(pc = 'PC9'), 
            tmp %>%
              filter(metric == 'tas_iasd') %>%
              select(region, metric, PC10) %>%
              rename(tas_iasd = PC10) %>%
              mutate(pc = 'PC10')
              ) ->
tas_iasd


tas %>%
  select(-metric) %>%
  left_join(pr %>% select(-metric), by =c('region', 'pc')) %>%
  left_join(tas_iasd %>% select(-metric), by =c('region', 'pc')) ->
  colored_pcs


colored_pcs$r <- scales::rescale(colored_pcs$tas_end_anomaly, to =c(0,255))
colored_pcs$g <- scales::rescale(colored_pcs$tas_iasd, to =c(0,255))
colored_pcs$b <- scales::rescale(colored_pcs$pr_end_anomaly_pct, to =c(0,255))
  
colored_pcs %>%
  left_join(as.data.frame(shp) %>% select(Acronym, mean_lon, mean_lat), by = c('region' = 'Acronym')) %>%
  rename(lon  = mean_lon, lat = mean_lat) %>%
  mutate(color = rgb(.$r, .$g, .$b, maxColorValue  = 255),
         color_order = as.integer(row.names(.))) ->
  pc_plotting

pc_plotting %>%
  select(color_order, color) %>%
  distinct() ->
  colors

shp %>% 
  left_join(pc_plotting, by = c('Acronym' = 'region'))->
  shp_pcs

p<- ggplot() +
  geom_sf(data = shp_pcs %>% filter(pc != 'PC10') %>% na.omit, aes(fill = as.factor(color_order)) ) +
  scale_fill_manual(values =colors$color)+
  facet_wrap(~pc, ncol=3) +
  theme(legend.position = 'none') 
p

#look at whether Precip inc or dec -> not just the PCs' precip entries >= 0
# because standardized/

as.data.frame(data_pca$center) %>%
  mutate(row_id = row.names(.)) %>%
  filter(grepl('pr_end_anomaly_pct', row_id)) %>%
  separate(row_id, into  = c('region', 'trash'), sep = '~') %>%
  select(-trash) ->
  pr_centers

as.data.frame(data_pca$scale) %>%
  mutate(row_id = row.names(.)) %>%
  filter(grepl('pr_end_anomaly_pct', row_id)) %>%
  separate(row_id, into  = c('region', 'trash'), sep = '~') %>%
  select(-trash) ->
  pr_scales




shp_pcs %>%
  left_join(pr_centers, by = c('Acronym'='region')) %>%
  left_join(pr_scales, by = c('Acronym'='region')) %>%
  mutate(precip_change = if_else( ((pr_end_anomaly_pct*`data_pca$scale`)- `data_pca$center`) >=0, 'inc', 'dec'))->
  shp_pcs2


ggplot() +
    geom_sf(data = shp_pcs2 %>% filter(pc != 'PC10') %>% na.omit, aes(#fill = as.factor(color_order),
                                             color = precip_change )) +
  #scale_fill_manual(values =colors$orig_color)+
  scale_color_manual(values = c('red', 'blue'))+
  facet_wrap(~pc, ncol=3) +
  theme(legend.position = 'none')